The most widely used techniques for geological surveying and hydrocarbon exploration are seismic methods. The seismic methods can image the structures of the sub-seafloor strata and reveal the location and shape of a potential reservoir, but face well documented difficulties in determining reservoir saturation. Conventionally, the solution to this is to drill a borehole into the reservoir. The costs of drilling a borehole offshore are expensive, often in tens of million dollars. Very recently, electromagnetic methods, for example controlled source electromagnetic methods (“CSEM”) and magnetotellurics (“MT”), have been developed for determining or mapping sub-seafloor resistivity variations. See, for example, U.S. Pat. No. 6,603,313 to Srnka. While seismic properties of hydrocarbon-filled and water-filled reservoirs do not differ significantly, their electromagnetic properties can be significantly different. For example, the resistivity difference between the two cases can be up to two orders in magnitude. Electromagnetic (“EM”) methods exploit these differences to predict the nature of a reservoir and save cost in hydrocarbon exploration. EM data inversion provides a technology to realize this exploitation in hydrocarbon exploration.
Geophysical data inversion is a procedure for obtaining earth models that satisfy measured data sets. The inversion process can provide physically meaningful information concerning both rock properties and earth structure, and therefore is a useful tool for the earth scientists. Inversion has been applied in global seismology, exploration seismic, potential field, and electromagnetic exploration. See, for example, Inversion of Geophysical Data, L. R. Lines, Ed., Society of Exploration Geophysicists (1988). 3D inversion of EM data can provide unique information related to reservoir location, shape and fluid properties. However, current 3D EM inversion schemes require expensive computer resources even to obtain low-resolution images.
The inversion process is closely related to forward modeling. Forward modeling uses a mathematical relationship (Maxwell's electromagnetic field equations for CSEM and MT) to simulate the earth's response for a given set of model parameters. Forward modeling can be written symbolically as d=F(m), where (for electromagnetic problems) m is a model of the earth's resistivity, F is known from Maxwell's equations for the EM fields, and d is a vector of response of the model m. Forward modeling provides a means to compute d for any model m. The inverse problem corresponding to this forward problem would be to find the set of all m that yield the given data d (a field or synthetic data set for inversion). It may be written, again symbolically, as m=F−1(d). This inverse operator F−1 is nonlinear, very complicated and non-unique for EM inversion. A simple and computationally tractable approach to the nonlinear multi-dimensional inverse problem is the linearized inversion. The nonlinear relationship between data and model in the forward problem is approximated by d=F(m0)+Gδm. The model update δm to a known (or guessed) model m0 can be obtained by solving a linear system Gδm=b, where G is the Jacobian matrix and b=d−F(m0) is the data residual. The model can be updated iteratively by adding δm to m0 until a satisfactory fit to the data has been obtained. The inverse problem and its solutions have been studied extensively (see, for example, R. L. Parker, Geophysical Inverse Theory (1994); W. Menke, Geophysical Data Analysis: Discrete Inverse Theory (1989); and A. Tarantola, Inverse Problem Theory, (1987)).
There are at least four major problems with EM inversions: The first problem is that many solutions are acceptable for inversion from a mathematical viewpoint (i.e., the non-uniqueness problem), especially when the data are limited and inaccurate. Mathematical approaches such as regularization are often implemented in inversion to mitigate some aspects of non-uniqueness. The second problem with multi-dimensional inversion is the cost. The linear equation system resulting from linearization is often very large for a multi-dimensional problem and requires the use of a supercomputer or massively paralleled computer system, particularly for CSEM data inversion. Newman and Alumbaugh inverted a synthetic data set of 12,600 source-receiver pairs with 1 frequency for a moderate 3D model with 29,971 cells. The processing time needed to produce a useful image was approximately 31 hours on the 1728-processor Intel Paragon, with 512 processors utilized (Newman and Alumbaugh, “Three-dimensional massively parallel electromagnetic inversion—I. Theory,” Geophys. J. Int. 128, 345-354 (1997)). It could take over a month to invert a large EM field data set for a large 3D model even on modern massively-parallel machines, which therefore limits 3D EM inversion application. The third problem with EM inversion is related to inversion resolution. Due to the diffusive nature of the EM field at low frequency, the resolution provided by EM is very low and cannot compete with seismic resolution. In general, a highly simplified picture (a blurred image) of 3D structures is all that can be obtained for EM inversion (West and Macnae, “Physics of Electromagnetic Induction Exploration Method,” in Electromagnetic Methods in Applied Geophysics (ed. M. N. Nambighian), Vol. 2, 5-45, Society of Exploration Geophysicists (1987)).
The fourth problem is related to model discretization in numerical modeling. Both multi-dimensional forward modeling and inversion are normally based on a discretized model. Discretization depends on the employed numerical modeling method (for example, the finite difference method, the integral equation method, the finite element method, or other method). Features and requirements on discretization for each method can easily be found in numerical modeling books. In principle, in order to model fine structures and achieve accurate forward modeling results, the model needs be discretized finely enough. FIG. 1 shows a uniform rectangular grid which is preferred for the finite difference method. This same fine discretization is typically used for both forward and inversion. That is, numerical methods of iteratively solving the inversion problem necessarily involve forward modeling at each iteration step, and it is typical to use the same discrete grid for each. The use of fine discretization in the inverse process has adverse effects: (1) it generates a huge system of linear equations (specially for the finite difference method), which needs a lot of computer resources for a reasonable turn-around time; (2) it may worsen the non-uniqueness because EM methods cannot resolve a small cell at depth due to lack of sensitivity (For example, cell j in FIG. 1). From another viewpoint, the contributions to response due to source TX measured at receiver RX from cell i at shallow depth and a deeper lying cell j of the same size are significantly different. In other words, receiver RX is much more sensitive to cell i than to cell j. Therefore, it is not optimal to treat them in the same way in inversion. The inversion non-uniqueness and the low EM resolution affect each other to make EM inversion much more challenging.
In reality, the unknown model m is a function of position, which is of infinite dimension, and the measurements d comprise only a finite collection of numbers with error, so that the inverse problem is not unique. As already mentioned non-uniqueness is a serious problem in inversion. A variety of approaches have been proposed to deal with the non-uniqueness problem. One approach has been to find localized averages that are shared by all models that are close enough to some reference model for a linearization approximation to hold (Backus and Gilbert, “The resolving power of gross earth data,” Geophy. J. R. astr. Soc. 16, 169-205 (1968); Parker, “The inverse problem of electromagnetic induction: existence and construction of solutions based on incomplete data,” J. Geophys. Res. 85, 4421-4428 (1970); Oldenburg, 1979, “One dimensional inversion of natural source magnetotelluric observations,” Geophysics 44, 1218-1244 (1979)). A second approach has been to find models minimizing some functional, particularly functionals that penalize roughness of the model (Tikhonov and Arsenin, Solutions of ill-posed problems, John Wiley and Sons (1977); Parker, “The theory of ideal bodies for gravity interpretation,” Geophy. J. R. astr. Soc. 42, 315-334 (1975); Constable et al., “Occam's inversion: a practical algorithm for generating smooth models from EM sounding data,” Geophysics 52, 289-300 (1987); and Smith and Booker, “Magnetotelluric inversion for minimum structure,” Geophysics 53, 1565-1576 (1988)). A third approach has been to assume prior knowledge of the distribution of likely models and find which of these models is most likely given a set of data (Franklin, “Well-posed stochastic extensions of ill-posed problems,” J. Math. Anal. Appl. 682-716 (1970); and Jordan and Franklin, “Optimal solutions to a linear inverse problem in geophysics,” Proc. Nat Acad. Sci. 68, 291-293 (1971)). A fourth approach (called joint inversion or cooperative inversion) has been to jointly invert various independent geophysical data sets (Vozoff and Jupp, “Joint inversion of geophysical data,” Geophy. J. R. astr. Soc. 42, 977-991 (1974); Savino et al., “Simultaneous inversion of multiple geophysical data sets for earth structure,” SEG 45th Annual International Meeting (1980); Lines et al., “Cooperative inversion of geophysical data,” Geophysics 53, 8-20 (1988); and Benech et al., “Joint inversion of EM and magnetic data for near-surface studies,” Geophysics 67, 1729-1739 (2002)). When dealing with specific or known structures derived from other geological and geophysical information, constrained inversion is employed to incorporate the specific structures into inversion. Wu incorporated structural constraints into model by inserting a discontinuous boundary within model; freezing the model at the specific nodes; and allowing different measures of model roughness in specified areas (“High resolution electromagnetic image of conductivity structure in the mid-lower crust and upper mantle—A magnetotelluric experiment conducted primarily in North Dakota,” Ph.D. Dissertation, Univ. of Washington (1994)). Structural information is also incorporated into seismic tomography (Grau and Lailly, “Sequential migration-aided reflection tomography: an approach to imaging complex structures,” Journal of Applied Geophysics 30; 75-87 (1993); Clapp, et al., “Incorporating geological information into reflection tomography,” Geophysics 69, 533-546 (2004)). The first three approaches use mathematic constraints to mitigate the non-uniqueness. Such constraints may not be consistent with the reality and therefore the inversion may provide an image inconsistent with the truth. The last two approaches use physical constraints from independent data sets which are consistent with the reality.
All of the preceding approaches result in a very large linear system to solve. Options to solve a very large linear system are very limited. A powerful computer is often needed in order to obtain results in a reasonable time. Nevertheless, a number of techniques have been developed to speed up computation at different stages of the inverse process. For example, more efficient optimization techniques such as non-linear conjugate gradient (NLCG) solver, multi-grid for modeling, approximate computation for sensitivity matrix, reciprocity application for source and receiver configurations, etc. All those techniques are helpful, but more improvements are needed to make 3D EM inversion a routine practice with reasonable demand on computer resources.
Non-uniqueness and resolution affect each other. Theoretically, the mathematical approaches to the non-uniqueness problem do not provide new information to enhance the resolution of a data set, rather than post constraints on model. However, they do affect the final image because of the implementation of mathematical constraints on the model. Joint inversion and constrained inversion not only mitigate the non-uniqueness but also enhance the resolution. Joint inversion results in a much larger linear system and therefore is more expensive to apply in practice. Even though constrained inversion utilizes some structural information in inversion, its inversion is still implemented to recover the geometry of structures as well as their physical properties. Therefore, the resolution provided by constrained EM inversion is largely limited by the diffusion nature of EM fields.
It is convenient and simple to use the same discretization in implementation of both forward modeling and inversion. No techniques that deal with the problem of using the fine discretization in inversion were found in EM inversion publications. In seismic tomographic inversion, the matrix transformation is used in order to use non-uniform grids, which is better for constructing a physical property model of a subsurface region (PCT International Publication No. WO2007/018869).
The present invention provides a faster method for inverting EM data with lower demand on computer resources for physically constrained solutions of high resolution.